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Neutron-rich matter at subnuclear densities may involve complex structures displaying a variety 
of shapes, such as spherical, slablike, and/or rodlike shapes. These phases of the nuclear pasta are 
expected to exist in the crust of neutron stars and in core-collapse supernovae. The dynamics of 
core-collapse supernovae is very sensitive to the interactions between neutrinos and nucleons/nuclei. 
Indeed, neutrino excitation of the low-energy modes of the pasta may allow for a significant energy 
transfer to the nuclear medium, thereby reviving the stalled supernovae shock. The linear response 
of the nuclear pasta to neutrinos is modeled via a simple semi-classical simulation. The transport 
mean-free path for ^ and r neutrinos (and antineutrinos) is expressed in terms of the static structure 
factor of the pasta, which is evaluated using Metropolis Monte Carlo simulations. 

PACS numbers: 



I. INTRODUCTION 

Neutron-rich matter may have a complex structure at 
densities just below that of normal nuclei. This is because 
all conventional matter is frustrated. While nucleons are 
correlated at short distances by attractive strong interac- 
tions, they are anti-correlated at large distances because 
of the Coulomb repulsion. Often these short and large 
distance scales are well separated, so nucleons bind into 
nuclei that are segregated in a crystal lattice. However, 
at densities of the order of 10'^'^— 10^** g/cm^ these length 
scales are comparable [J- Competition among these in- 
teractions {i.e., frustration) becomes responsible for the 
development of complex structures with many possible 
nuclear shapes, such as spheres, cylinders, flat plates, as 
well as spherical and cylindrical voids . The term pasta 
phases has been coined to describe these complex struc- 
tures T|, and many calculations of their ground-state 
structure have already been reported 0, |2j |3j IS - 
While the study of these pasta phases is interesting in its 
own right, it becomes even more so due to its relevance 
to the structure of the inner crust of neutron stars and 
to the dynamics of core-collapse supernovae. 

Frustration, a phenomenon characterized by the exis- 
tence of a very large number of low-energy configurations, 
emerges from the impossibility to simultaneously mini- 
mize all "elementary" interactions. Should a proton in 
the pasta join a nuclear cluster to benefit from the nuclear 
attraction or should it remain well separated to minimize 
the Coulomb repulsion? Frustration, a term that appears 
to have been coined in the late seventies 0, Q , is preva- 
lent in complex systems ranging from magnetism ]^ [To| 
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to protein folding In condensed-matter systems, 

frustration dates back to the 1950 study of Ising an- 
tiferromagnets on triangular lattices by Wannier |l2j |. 
Three antiferromagnetically coupled spins fixed to the 
sites of an equilateral triangle cannot minimize all inter- 
actions simultaneously: once two spins are anti-aligned, 
the third one cannot be antiparallel to both of them. Fur- 
ther, in Ref. it has been shown that finding the true 
ground state — among the many metastable states — of a 
spin glass shares features in common with NP- complete 
problems, such as the traveling salesman problem of fame 
in the theory of combinatorial optimization . Finally, 
because of the preponderance of low-energy states, frus- 
trated systems display unusual low-energy dynamics. 

Because of the complexity of the system, almost no 
work has been done on the low-energy dynamics of the 
pasta or on its response to weakly-interacting probes. In 
this paper we study the excitations of the pasta via a 
simple semi-classical simulation similar to those used to 
describe heavy-ion collisions. Heavy-ion collisions can 
produce hot, dense matter. However, by carefully heat- 
ing the system and allowing it to expand, heavy-ion col- 
lisions can also study matter at low densities. Multi- 
fragmentation, the breakup of a heavy ion into several 
large fragments, shares many features with pasta for- 
mation, as they are both driven by the same volume, 
surface, and Coulomb energies. There have been sev- 
eral classical |l5l Il6j and quantum-molecular-dynamics 
(QMD) simulations of heavy-ion collisions. These 
same approaches may be applied to the nuclear pasta by 
employing a simulation volume and periodic boundary 
conditions. One great advantage of such simulations is 
that one can study pasta formation in an unbiased way 
without having to assume particular shapes or configura- 
tions from the outset. While QMD has been used before 
to study the structure of the pasta 18, 19, 20], no calcu- 
lations of its linear response to weakly-interacting probes 
{e.g., neutrinos) have been reported. 
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Neutrino interactions are crucial in the dynamics of 
core-collapse supernovae because neutrinos carry 99% of 
the energy. Neutrinos are initially trapped due to the 
large coherent neutrino-nucleus elastic scattering cross 
section. This trapping is important for the electron-per- 
baryon fraction Ye of the supernova core, as it hinders 
any further conversion of electrons into neutrinos. How- 
ever, with increasing density Coulomb interactions be- 
tween ions lead to significant ion screening of neutrino- 
nucleus elastic scattering (21| . As the density is increased 
still further, the ions react to form nuclear pasta and one 
needs to calculate neutrino-pasta interactions. 

In the pasta phase one can have coherent neutrino scat- 
tering from density contrasts, such as Swiss cheese like 
voids. Critical fluctuations could significantly increase 
the cross section, thereby greatly reducing the neutrino 
mean-free path. As the existence of many low-energy 
configurations is the benchmark of frustrated systems, 
one expects large configuration mixing among the vari- 
ous different pasta shapes. This mixing often leads to 
interesting low-energy collective excitations. As it will 
be shown later (see Fig. ^ at subnuclear densities of 
the order of 10^^ g/cm'^, the pasta resembles a collec- 
tion of spherical neutron-rich nuclei embedded in a di- 
lute neutron gas. Neutron-rich nuclei with large neu- 
tron skins have Pygmy giant resonances, involving col- 
lective oscillations of the neutron skin against the sym- 
metric core [23, li^l- We expect that the soft neutron- 
rich pasta will have many low-energy collective oscil- 
lations. This could provide important physics that is 
presently missing from core-collapse supernovae simula- 
tions. Neutrino excitation of the low-energy pasta modes 
may allow for a significant energy transfer to the nu- 
clear medium, potentially reviving the stalled supernovae 
shock. To our knowledge, there have been no calculations 
of these effects. Note, however, that Reddy, Bertsch, and 
Prakash j24| have found that coherent neutrino scatter- 
ing from a nonuniform kaon condensed phase greatly de- 
creases the neutrino mean free path. 

Present models of the equation of state for supernovae 
simulations, such as that of Lattimer and Swesty [25l| . 
describe the system as a liquid drop for a single repre- 
sentative heavy nucleus surrounded by free alpha parti- 
cles, protons, and neutrons. One then calculates neu- 
trino scattering from these constituents — by arbitrarily 
matching to a high-density uniform phase [2g . Unfortu- 
nately, this approximation is uncontrolled as it neglects 
many important interactions between nuclei. By simulat- 
ing the pasta phase directly in the nuclcon coordinates, 
one hopes to improve on this matching and to understand 
its limitations. 

There is a duality between microscopic descriptions of 
the system in terms of nucleon coordinates and "macro- 
scopic" descriptions in terms of effective nuclear degrees 
of freedom. Thus, a relevant question to pose is as fol- 
lows: when does a neutrino scatter from a nucleus and 
when does it scatter from an individual nucleon? At the 
Jefferson Laboratory a similar question is studied; when. 



i.e., at what momentum transfer, does a photon couple 
to a full hadron and when to an individual quark? Mod- 
els of the quark/hadron duality have provided insight on 
how descriptions in terms of hadron degrees of freedom 
can be equivalent to descriptions in terms of quark coor- 
dinates |27j |. Here we are interested in nucleon/nuclear 
duality, that is, how can nuclear models incorporate the 
main features of microscopic nucleon descriptions? 

The manuscript has been organized as follows. In 
Sec.mthe semi-classical formalism is introduced. A very 
simple (perhaps minimal) model is employed that con- 
tains the essential physics of frustration. The linear re- 
sponse of the pasta to neutrino scattering, in the form of 
a static structure factor, is discussed in Sec. IIIII Results 
are presented in Sec. IIVI while conclusions and future 
directions are reserved to Sec. |3 



II. FORMALISM 

In this section we introduce a classical model that while 
simple, contains the essential physics of frustration. That 
is, competing interactions consisting of a short-range nu- 
clear attraction and a long-range Coulomb repulsion. We 
model a charge-neutral system of electrons, protons, and 
neutrons. The electrons are assumed to be a degenerate 
free Fermi gas of density pe = Pp and the nucleons inter- 
act via a classical potential. The only quantum aspects of 
the calculation are the use of an effective temperature and 
effective interactions to simulate effects associated with 
quantum zero-point motion. Of course more elaborate 
models are possible and these will be presented in future 
contributions. For these first simulations we adopted a 
very simple version that displays the essential physics 
of nucleons clustering into pasta in a transparent form. 
Moreover, this simple model facilitates simulations with 
a relatively large numbers of particles, a feature that is 
essential to estimate and control finite-size effects. 

The total potential Vtot energy is assumed to be a sum 
over two-body interactions Vij of the following form: 



(1) 



where the "elementary" two-body interaction is given by 

V{i,j) = ae--"^^^+[b + cT,{{)T,{j)\e'^"^/^'' + V,{i,j) . 

(2) 

Here the distance between the particles is denoted by 
Tij — Ir— fj l and the isospin of the j'th particle is T2(j) = 1 
for a proton and t^(j) = —1 for a neutron. The model 
parameters a, 6, c, and A will be discussed below. Suf- 
fices to say that the above interaction includes the char- 
acteristic intermediate-range attraction and short-range 
repulsion of the nucleon-nucleon {NN) force. Further, 
the isospin dependence of the potential insures that while 
pure neutron matter is unbound, symmetric nuclear mat- 
ter is bound appropriately. Finally, a screened Coulomb 



interaction of the following form is included: 

=2 



function as follows: 



(3) 



where Tp{j) ^ {l + Tz{j))/2 and A is the screening length 
that results from the slight polarization of the electron 
gas. That is, the rclativistic Thomas-Fermi screening 
length is given by 



A = - ( kpJkl + ml 



-1/2 



(4) 



Note that the electron Fermi momentum has been defined 
by fcp = (Stt^Pe)-'^/^ and rrie is the electron mass |28ll29| . 
Unfortunately, while the screening length A defined above 
is smaller than the length L of our simulation box, it 
is not significantly smaller (unless a prohibitively large 
number of particles is used). Therefore, to control finite- 
size effects we were forced to arbitrarily decrease the 
value of A (see Sec. lIV|l . 

The simulations are carried out in a canonical ensem- 
ble with a fixed number of particles A at a temperature 
T . The volume y at a fixed baryon density p is simply 
V = A/ p. To minimize finite-size effects we use periodic 
boundary conditions, so that the distance is calcu- 
lated from the a;, y, and z coordinates of the ith and jth 
particles as follows: 



(5) 



where the periodic distance, for a cubic box of side L — 
is given by 



[/]=Min(|/|,L-|/|). 



(6) 



The potential energy defined in Eq. is indepen- 
dent of momentum. Therefore, the partition function for 
the system factors into a product of a partition function 
in momentum space — that plays no role in the compu- 
tation of momentum-independent observables — times a 
coordinate-space partition function of the form 



Z{A, T, y) = j (fri ■ ■ ■ (frA exp (- 



■Vtot/T) . (7) 



Note that the 3-diniensional integrals are performed over 
the simulation volume V . 

The average energy of the system {E) — {K) + {Vtot) is 
made of kinetic {K) and potential (T4ot) energy contri- 
butions. As the (momentum-independent) interactions 
have no impact on momentum-dependent quantities, the 
expectation value of the kinetic energy reduces to its clas- 
sical value, that is. 



(K) = -AT 



(8) 



In turn, the expectation value of the potential energy 
may be computed from the coordinate-space partition 



{Vtot) = 



1 



Z{A,T,V) 



dPri ■ ■ ■ dpTAVtot exp (-14ot/T) 



(9) 

In summary, a classical system has been constructed 
with a total potential energy given as a sum of two- 
body momentum- independent interactions [see Eq. Q]. 
Any momentum-independent observable of interest can 
be calculated from the partition function [Eq. (|7I)], which 
we evaluate via Metropolis Monte-Carlo integration |30l| . 
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FIG. 1: Energy per particle for symmetric (dashed line) and 
for pure-neutron matter (solid) versus baryon density nh at a 
temperature of r = l MeV. 



We now return to discuss the choice of model param- 
eters. The constants a, 6, c, and A in the two-body in- 
teraction Eq. 121 were adjusted — approximately — to re- 
produce the following bulk properties: a) the saturation 
density and binding energy per nucleon of symmetric nu- 
clear matter, b) (a reasonable value for) the binding en- 
ergy per nucleon of neutron matter at saturation density, 
and c) (approximate values for the) binding energy of a 
few selected finite nuclei. The temperature was arbitrar- 
ily fixed at 1 MeV for all the calculations. Note that 
the parameter set employed in these calculations (and 
displayed in Table has yet to be carefully optimized. 
We reiterate that for these first set of simulations, the 
interaction is sufficiently accurate to describe the essen- 
tial physics of the pasta. Indeed, this is illustrated in 
Fig. n and Table ITU In Fig. ^ the average potential en- 
ergy versus density at a temperature of T = 1 MeV is 
displayed for a simulation of symmetric nuclear matter 
containing A = 400 particles and, as is customary, assum- 
ing no Coulomb interactions. Also shown in the figure is 
the potential energy for pure neutron matter calculated 
with N — 200 particles. In the case of finite nuclei (also 
calculated at T = 1 MeV) the full Coulomb interaction is 
included using a screening length A much larger than the 
resulting root-mean-square radius of the nucleus. Simula- 
tions based on a Metropolis Monte-Carlo algorithm were 
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TABLE I: Model parameters used in the calculations. 



a 


b 


c A 


110 MeV 


-26 MeV 


24 MeV 1.25 fm^ 


TABLE IL Binding energy per nucleon for various closed shell 
nuclei. Note that all energies are in MeV and that the Monte- 
Carlo results include only (Vtot). 


Nucleus 


Monte Carlo (Vtot) 


Experiment 


16q 


-7.56 ±0.01 


-7.98 


*«Ca 


-8.75 ±0.03 


-8.45 


90Zr 


-9.13 ±0.03 


-8.66 


208pb 


-8.2 ±0.1 


-8.45 



used to compute the average potential energy, starting 
with nucleons distributed uniformly in a sphere with a 
radius comparable to the expected size of the nucleus; 
this sphere was placed in the center of a very large box. 
Results of the simulations and comparison with experi- 
mental values have been collected in Tabled] Note that 
the simulation results are for the potential energy only. If 
the kinetic energy per nucleon {3T/2) is added to these 
values, the nuclei in Table ^ would be slightly under- 
bound. Furthermore, finite nuclei are only metastable 
in this semiclassical approximation. Nucleons can evap- 
orate over a very long time scale. However, this is not 
expected to significantly impact the pasta phases since 
these already have free nucleons. 

III. NEUTRINO SCATTERING 

The model is used to describe neutrino scattering from 
nonuniform neutron-rich matter. As neutrino interac- 
tions play a predominant role in core-collapse supernovae, 
one is interested in understanding how the neutrinos dif- 
fuse and how do they exchange energy. In this first paper 
we focus on the transport mean free path for and Vr, 
which lack charged-current interactions at low energies. 
Their mean-free path is dominated by neutral current 
neutrino-nucleon scattering. 

The free-space cross section for neutrino-nucleon elas- 
tic scattering is given by 

^-^K(3-cos0)±c^(l±cos^)] , (10) 

where Gp is the Fermi coupling constant, E^, the neu- 
trino energy and 9 the scattering angle. Note that this 
equation neglects weak magnetism and other corrections 
of order E^/M, with M the nucleon mass |3]| . 

In the absence of weak magnetism, the weak neutral 
current of a nucleon has axial- vector (75 7^^) and vector 
7;i contributions. That is, 

= Ca 757m + Ci,7m ■ (11) 



The axial coupling constant is, 

Ca=±|^ (5a = 1.26), (12) 

with the + sign for neutrino-proton and the — sign for 
neutrino-neutron scattering. The weak charge of a proton 
Ct, is suppressed by the weak-mixing (or Weinberg) angle 
sin^ 6lw = 0.231, 

c„ = i - 2 sin^ 0w = 0.038 w . (13) 

In contrast, the weak charge of a neutron is both large 
and insensitive to the weak-mixing angle: c^ = — 1/2. The 
transport mean-free path At is inversely proportional to 
the transport cross section at and is given by the follow- 
ing expression: 

= / dn^il - cose) = (5c^ + c^) . (14) 

The weighting factor (1— cos0) included in the definition 
of the transport cross section favors large-angle scatter- 
ing, as momentum is transferred more efficiently into the 
medium. As a result, the axial-vector contribution 
dominates the cross section. Assuming that the scat- 
tering in the medium is the same as in free space, the 
transport mean-free path becomes 

At = (pX + Pnar)-' . (15) 

Here Pp{pn) is the proton (neutron) density and aficrt) 
is the transport cross section for scattering from a pro- 
ton(neutron). 

If nucleons cluster tightly into nuclei or into pasta, 
then the scattering from different nucleons could be co- 
herent. This will significantly enhance the cross section 
as it would be proportional to the square of the number of 
nucleons [s^. The contribution from the vector current 
is expected to be coherent. Instead, the strong spin and 
isospin dependence of the axial current should reduce its 
coherence. This is because in nuclei — and presumably in 
the pasta — most nucleons pair off into spin singlet states 
(note that in the nonrelativistic limit the nucleon axial- 
vector current becomes 757T2 — ctt^). Therefore, in 
this paper we focus exclusively on coherence effects for 
the vector current. 

Coherence is important in x-ray scattering from crys- 
tals. Because the x-ray wavelength is comparable to the 
inter-particle spacing, one needs to calculate the rela- 
tive phase for scattering from different atomic planes and 
then sum over all planes. Neutrino-pasta scattering in- 
volves a similar sum because the neutrino wavelength is 
comparable to the inter-particle spacing and even to the 
inter-cluster spacing. Therefore, one must calculate the 
relative phase for neutrino scattering from different nu- 
cleons and then add their contribution coherently. This 
procedure is embodied in the static structure factor S{q). 
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The dynamic response of the system to a probe of mo- 
mentum transfer q and energy transfer w > that couples 
to the weak charge density p(q) is given by 



5(q,c.) 



^1 



(16) 



where uj„ is the energy difference between the excited 
state l^'n) and the ground state |4'o)- In linear response 
theory, namely, assuming that the process can be treated 
in lowest order (an excellent approximation for neutrino 
scattering) the cross section can be directly related to 
the dynamic response of the system. In the case that 
the individual excited states may not be resolved, then 
one integrates over the energy transfer to to obtain the 
static structure factor. Here we define the static structure 
factor per neutron as follows: 



1 r°o 1 
^(q) - -TT / Siq,^)dLU = - ^ |(*„|p(q)|*o 



with the weak vector-charge density given by 



N 



P(q) = exp(«q • r,) 



(17) 



(18) 



where the sum in Eq. (|18|l is only over neutrons. 

The cross section per neutron for neutrino scattering 
from the whole system is now given by 



cos 6) 



(19) 



Note that the weak charge of the nucleon {cy « for pro- 
tons and Cy = —1/2 for neutrons) has been incorporated 
into the above cross section, so that the normalization 
of the weak vector-charge density is p{q = 0) = N. Fur- 
ther, Eq. H19|l is the cross section per neutron obtained 
from Eq. lO (with Ca = 0) multiplied by S(q). This 
indicates that S(q) embodies the effects from coherence. 
Finally, note that the momentum transfer is related to 
the scattering angle through the following equation: 



g2 2£;^(i _cos6i) 



(20) 



Two assumptions have been made in the derivation of 
Eq. (|19|l . First, no contribution from the axial current 
to the cross section has been included, because nucleons 
pair into spin-zero states. Second, the excitation energy 
transferred to the nucleons is small and we have summed 
over all possible excitation energies. 

The static structure factor has important limits. A 
particularly useful form in which to discuss them invokes 
completeness on Eq. (|17ll . That is. 



:*o|/5t(q)/5(q)|vI/o) - (*o|/5(q)|* 



(21) 



The last term in the above expression represents the elas- 
tic form factor of the system, which only contributes at 
q = 0. In the limit oi q ^ 0, the weak charge den- 
sity [Eq. H18|l ] becomes the number operator for neutrons 
p{q — 0) — N, so that the static structure factor reduces 
to, 



5(<Z = 0) = l((7V2)-(iV)2 



(22) 



Thus, the g— >0 limit of the static structure factor is re- 
lated to the fluctuations in the number of particles, or 
equivalently, to the density fluctuations. These fluctua- 
tions are, themselves, related to the compressibility and 
diverge at the critical point js^l- To discuss the large q 
limit, Eq. (fTH|) is substituted into Eq. lf?T)) to yield 



Siq)-- 




^'o|exp(iq- rjj)|*o)- (*o|/5(q)|*o 



(23) 

In the 5 — > oo limit, all the terms in the sum with ij^j, as 
well as the second term in the above expression, oscillate 
to zero. This only leaves the i = j terms, which there are 
N of them so that 



Siq 



oo 



1 



(24) 



This result indicates that if the neutrino wavelength is 
much shorter than the inter-particle separation, the neu- 
trino only resolves one nucleon at a time. This corre- 
sponds to quasielastic scattering where the cross section 
per nucleon in the medium is the same as in free space. 

One can calculate the static structure factor from the 
neutron-neutron correlation function which is defined as 
follows: 



1 



N 

— F< 

Npn ^ 



(25) 



The two-neutron correlation function "asks" (and "an- 
swers" ) the following question: if one "sits" on a neutron, 
what is the probability of finding another one a distance 
|r| away. The correlation function is normalized to one 
at large distances g{r — oo) — 1; this corresponds to 
the average density of the medium. The static structure 
factor is obtained from the Fourier transform of the two- 
neutron correlation function. Comparing with Eq. I|23l) 
this yields. 



Siq) = l + pnj rfV(5(r) - l) exp(zq • r) 



(26) 



The i—j terms in Eq. (|23|l gives the leading 1 in the above 
expression, while the elastic form factor |(^'o|/5(q)|'I'o)P 
yields the —1 in the integrand of Eq. I|26|) . 

To obtain the transport cross section we proceed, as in 
Eq. H1(J|) , to integrate the angular- weighted cross section 
da/di^{l — cos 6) over all angles. That is. 



^t-^/ dn^a- cose) ^{S{E,))a: 



(27) 
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Note that the free neutron cross CF^ follows directly from 
Eq. UH) in the limit of Ca = 0, 



Further, the angle averaged static structure factor has 
been defined as follows: 

{S{E,)) = ^ y dx{\ - x^)S{q{x, E,)) , (29) 

where the static structure factor S{q) depends on neu- 
trino energy and angle through Eq. it^ . that is, = 
2El{l — x) with x = cos9. Using Eq. and switching 
the orders of integration, this can be written as 

(SiE,)) = 1 + ^ drf{2E,r)[g{r) - l) . (30) 

Note that the following function has been introduced 

f{t) = 72(cost + tsint -!)/<'' -6(5 cost + < sin i+l)/i^ . 

(31) 

Finally, the transport mean free path (Af = l/<7tp„) is 
given by 

Ar^=a,V(5(i?.)> . (32) 

In this way {S{Ey)) describes how coherence modifies the 
mean-free path. In the next section, simulation results for 
the two-neutron correlation function, the static structure 
factor, and the angle averaged static structure factor will 
be presented. 



This value for A is adopted hereafter for all of our sim- 
ulations (see also This smaller screening length 
decreases slightly the Coulomb interaction at large dis- 
tances, which could promote the growth of slightly larger 
clusters. However, we do not expect this decrease in A to 
qualitatively change our results. We note that A is still 
larger than the typical size of a heavy nucleus. 




I \ \ \ \ 1 fm 

15 30 45 60 75 

FIG. 2: (Color online) Monte-Carlo snapshot of a configura- 
tion of A'' = 3200 neutrons (light gray circles) and Z = 800 pro- 
tons (dark red circles) at a baryon density of p = 0.01 fm~^, 
a temperature of T = 1 MeV, and an electron fraction of 
Fe = 0.2. 3D imaging courtesy of the FSU Visualization Lab- 
oratory. 



IV. RESULTS 

In this section we present our simulation results. As an 
example of a typical low-density condition we consider a 
subnuclear density of p = 0.01 fm^"^ (about 1/16 of nor- 
mal nuclear density), a temperature T=l MeV, and an 
electron fraction of Ye = 0.2. Proto-neutron stars have 
electron fractions that start out near 1/2 and drop with 
time, so Ye = 0.2 represents a typical neutron-rich condi- 
tion. Monte-Carlo simulations for a total of A = 4000 par- 
ticles (TV = 3200, Z = 800) have been performed. Because 
of the many competing minima, a significantly larger sys- 
tem would take an unreasonably long time to thermalize 
on a modest work station (see details below). 

The simulation volume for the above conditions con- 
sists of a cube of length L — 73.7 fm. While this value is 
larger than the electron screening length A = 26.6 fm [see 
Eq. Q], it is not sufficiently larger. Indeed, to minimize 
finite-size effects in a simulation with periodic boundary 
conditions one would like exp(— L/(2A)) ^ 1. Clearly, 
this condition is not adequately satisfied. Therefore, 
in an effort to minimize the contamination from finite- 
size effects, we reduce the electron screening length — 
arbitrarily — to the following value: 

A = 10 fm . (33) 



By far, the most time consuming part of the simula- 
tion is producing suitable initial conditions. The simu- 
lations are started with the A = 4000 nucleons randomly 
distributed throughout the simulation volume. Next, we 
perform a total of about 325,000 Metropolis sweeps start- 
ing at the higher temperature of T = 2 MeV and reduc- 
ing the temperature until eventually reaching the target 
temperature of T = 1 MeV, in a "poor's-man" attempt 
at simulated annealing. Note that a Metropolis sweep 
consists of a single trial move for each of the A = 4000 
particles in this system. We call this procedure cooking 
the pasta. 

Results in this section are based on a statistical av- 
erage of the final 50,000 sweeps. This yields a potential 
energy of — 5.385±0.003 MeV/A. A sample configuration 
of the 4000 particles is shown in Fig. 13 The protons are 
strongly correlated into clusters ( "nuclei" ) as are a large 
number of neutrons. In addition, there is a low density 
neutron gas between the clusters. At this density it may 
be reasonable to think of the system as a high-density 
liquid of "conventional" nuclei immersed in a dilute neu- 
tron gas. A great virtue of the simulation is that one does 
not have to arbitrarily decide which nucleons cluster in 
nuclei and which ones remain in the gas. These "deci- 
sions" are being answered dynamically. Further, one can 
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calculate modifications to nuclear properties due to the 
interactions. In a future work we plan to compare our 
simulation results to some conventional nuclear models. 

Protons moving between clusters face a Coulomb bar- 
rier. This may inhibit the thermalization process and 
with it the formation of larger clusters. This could in- 
crease our results for S{q). To test the thermalization of 
the pasta, our Metropolis Monte Carlo configuration was 
evolved further via molecular dynamics for a total time 
of 46,500 fm/c. The molecular dynamics calculations will 
be described in future work. This led to an increase in the 
peak of S{q) by only about 10 percent. Although these 
molecular dynamics results did not reveal a large secu- 
lar change in the system with time, we caution that our 
cooking procedure may not have converged to the true 
thermal-equilibrium state. The static structure factor 
S{q) may still change with additional Metropolis Monte 
Carlo or molecular dynamics evolution. 




1 L-l I I I \ I \ I L 

10 20 30 40 



r(fm) 

FIG. 3: Neutron-neutron correlation function at a tempera- 
ture of r = 1 MeV, an electron fraction of Ye —0.2, and a 
baryon density of p = 0.01 fm~'^. These results show large 
finite-size effects beyond r — L/2 — 36.9 fm. 

The neutron- neutron correlation function g(r) is dis- 
played in Fig. 121 The two-neutron correlation function is 
calculated by histograming the relative distances between 
neutrons. The correlation function is very small at short 
distances because of the hard core in our NN interac- 
tion. At intermediate distances g{r) shows a large broad 
peak between r = 2 fm and r ~ 10 fm. This corresponds 
to the other neutrons bound into a cluster. Superim- 
posed on this broad peak one observes three (or four) 
sharper peaks corresponding to nearest, next-to- nearest, 
and next-to-next-to- nearest neighbors. These structures 
describe two-neutron correlations within the same clus- 
ter. At larger distances, between 10 and 20 fm, the cor- 
relation function shows a modest dip below one, suggest- 
ing that the attractive NN interaction has shifted some 
neutrons from larger to smaller distances in order to form 



the clusters. Alternatively, Coulomb repulsion makes it 
less likely to find two clusters separated by these radii. 
Finally, there is an abrupt drop in g(r) at a distance 
corresponding to half the value of the simulation length 
(r = L/2 = 36.9 fm) caused by finite-size effects. We note 
that under our assumptions of periodic boundary condi- 
tions, g{r) is identically zero for r>V^L/2. 




FIG. 4: (Color online) Monte-Carlo snapshot of a configura- 
tion oi N = 3200 neutrons (light gray circles) and Z — 800 pro- 
tons (dark red circles) at a baryon density of p = 0.025 fm~^, 
a temperature of T = 1 MeV, and an electron fraction of 
Ye — 0.2. 3D imaging courtesy of the FSU Visualization Lab- 
oratory. 

Increasing the density can change the nature of the 
pasta. Figure 0] shows a sample configuration of 4000 
particles at a density of 0.025 fm~'^. Note that although 
the density has increased, both the temperature and the 
electron fraction have remained fixed at T = 1 MeV and 
Ye — 0.2, respectively. At a density of p = 0.025 fm"'^ 
(about 1/6 of normal nuclear-matter saturation density) 
the spherical clusters of Fig. |21 start to coalesce into 
cylindrical-like structures. The system appears to be 
segregated into high-density regions of cylindrical nuclei 
immersed in a dilute neutron gas. As we continue to 
perform additional simulations, high-quality renderings 
of nucleon configurations for a variety of densities, tem- 
peratures, and electron fractions will be developed. 

To compute the static structure factor S{q), one nu- 
merically transforms the two-neutron correlation func- 
tion, as indicated in Eq. H26|) . However, because of finite- 
size effects we truncate the integral at r = L/2 and assume 
g{r) = l for r> L/2. This procedure yields the results dis- 
played in Fig. [S] There is a modest peak in S{q — 0) due 
to density fluctuations. Of course, the number of neu- 
trons in our simulation remains fixed, yet fluctuations 
can take neutrons across the r — L/2 cutoff and these 
fluctuations will contribute to the value of S{q) at ^ = 0. 

The error bars in Fig. [51 are statistical only, based on 
the last 50,000 sweeps. We caution that there may be 
finite-size effects at small momentum transfers. The box 
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FIG. 5: Static structure factor S{q) versus momentum trans- 
fer q at a temperature of r = 1 MeV, an electron fraction of 
Ye = 0.2, and a baryon density of p = 0.01 fm~^. The error 
bars are statistical only. Finite size effects may be important 
at small q as indicated by the dotted error bars. 



size for our simulation at p = 0.01 fm ^ is iy = 73.7 fm. 
This corresponds to a minimum momentum transfer of 

<Zmin ~ Y " ^-^^^ ^™ ' • ^^^^ 

Momentum transfers smaller than Qmin correspond to 
wavelengths larger than the simulation volume, so our 
results for q< 2gmin may be sensitive to finite-size effects. 
Indeed, for q < 2q^i^ the static structure factor was ob- 
served to change significantly from one Metropolis run 
to the next. To indicate the sensitivity of our results to 
finite-size effects, Fig. displays the static structure fac- 
tor in the q< 2qmin region with dotted error bars. Note 
that the point q — gmin has been signaled out in Fig. [51 
to indicate that it is more stable from one Metropolis 
run to the next, because the weak vector-charge density 
[Eq. H18|l ] evaluated at q = qminq is invariant under a 
translation of the system by a distance L along q. 

The static structure factor displays a large peak at 
qwO.3 fm^^, corresponding to coherent scattering from 
many neutrons bound into a single cluster. At smaller 
momentum transfers, g « 0.2 fm~^, S{q) decreases be- 
cause of ion screening. Here the neutrino wavelength is 
so long that it probes multiple clusters. These other clus- 
ters screen the weak charge and reduce the response. At 
momentum transfers larger than gwO.B fm~^, the static 
structure factor decreases with increasing q. This is the 
effect of the cluster form factor. As the momentum trans- 
fer increases the neutrino can no longer scatter coherently 
from all the neutrons in a cluster because of the clusters 
extended size. Thus, the observed peak in S{q) develops 
as a trade off between ion screening, which favors large 
q, and the cluster form factor, which favors small q. 



In summary, one can divide the response of the pasta 
into the following regions. At low- momentum transfers 
{q<0.2 fm"^) the response is dominated by ion screen- 
ing and density fluctuations. For momentum transfers 
in the region q = 0.2 — 0.4 fm^^ one observes coherent 
scattering from the pasta. At the larger momenta of 
q = 0.4—1 fm~^, the falling response reflects the pasta 
form factor. Finally, the large momentum transfer region 
above q=l fm~^ corresponds to quasielastic scattering 
from nearly free neutrons, as S{q — > oo) = 1 [see Eq. (|24(l ]. 




"o 20 40 60 80 100 

(MeV) 

FIG. 6: Angle averaged static structure factor {S{E^)) ver- 
sus neutrino energy E^, at a temperature of T = 1 MeV, 
an electron fraction of Ye = 0.2, and a baryon density of 
p = 0.01 fm-^ The error bars are statistical only, see text. 



The angle averaged structure factor {S{E^)) [defined 
in Eq. (|29() ] is shown in Fig. |B| Note that the integral 
in Eq. (|30|l was also truncated at r = L/2 because of 
finite-size effects. The averaged structure factor shows 
a broad peak for neutrino energies near 40 MeV. In- 
deed, the transport cross section is significantly enhanced 
by coherence effects for neutrino energies from about 20 
to 80 MeV. The impact of this coherence on the neu- 
trino mean-free path [Eq. H32|) ] is displayed in Fig. [7| 
Also shown in the figure is the mean-free path obtained 
by ignoring coherence effects by setting {S{E^)) = 1 in 
Eq. (|32|) . Coherence significantly reduces the mean- free 
path for neutrino energies in the range E^, = 15—120 MeV. 
Again, finite size effects may be important for low neu- 
trino energies. 



V. CONCLUSIONS 

Neutron-rich matter is expected to have a complex 
structure at subnuclear densities. Complex pasta phases 
may result from frustration through the competition be- 
tween an attractive nuclear interaction and the Coulomb 
repulsion. Neutrino interactions with the pasta may 
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FIG. 7: Transport mean-free path At for i/^ or Ut versus neu- 
trino energy Ei, at a baryon density of p = 0.01 fm~'^, a tem- 
perature of T = 1 MeV, and an electron fraction of Ye =0.2. 
The solid line and error bars include full coherence effects 
while the dashed line is obtained by using {S{E)} = 1 in 
Eq. 1321 . Error bars are statistical only. Finite size effects 
may be important for low as indicated by the dotted error 
bars. 



be important for properties of core-collapse supernovae, 
such as its electron fraction. 

In this work we have employed a semi-classical model 
to simulate the dynamics of the pasta phase of neutron- 
rich matter. Although our model is very simple, it 
nonetheless retains the crucial physics of frustration. Us- 
ing a Metropolis Monte-Carlo algorithm, the partition 
function was computed for a system of 4000 nucleons at 
a given temperature and density. We find that almost all 
protons and most of the neutrons cluster into "nuclei" 
that are surrounded by a dilute neutron gas. 

Observables computed in our simulations included the 
neutron-neutron correlation function. This calculation 
was implemented by constructing a histogram of all rel- 
ative neutron distances. The two-neutron correlation 
function g{r) gives the probability of finding a neutron 
at a distance r away from a reference neutron. A large 
peak in g{r) at intermediate distances (r = 2 — 10 fm) is 
found. This reflects the presence of the other neutrons 
in the cluster. 

The static structure factor S{q), a fundamental ob- 
servable obtained from the Fourier transform of the two- 
neutron correlation function, describes the degree of co- 
herence for neutrino-nucleon elastic scattering. For small 
momentum transfers, S'(q) describes density fluctuations 
and ion screening. In this region the neutrino wave- 
length is longer than the average inter-cluster separation, 
thereby allowing other clusters to screen the weak charge 
of a given cluster. At momentum transfers of approx- 
imately q = 0.2 — 0.4 fm~^, the static structure factor 
develops a large peak, associated to the coherent scatter- 



ing from all the neutrons in the cluster. This coherence 
is responsible for a significant reduction in the neutrino 
mean-free path. To our knowledge, these represent the 
first consistent calculation of the neutrino mean-free path 
in nonuniform neutron-rich matter. 

However, much remains to be done. First, one needs to 
focus on the thermalization of our simulations. It is diffi- 
cult to insure that the system has reached thermal equi- 
librium because the Coulomb barrier hinders the motion 
of individual protons. Second, one must further investi- 
gate the impact of finite-size effects and the simple treat- 
ment of long-range Coulomb interactions on our simula- 
tions. This may require simulations with larger numbers 
of particles, as it is difficult to fit a long-wavelength neu- 
trino into the present simulation volume. Third, while 
we have focused here on the vector part of the weak- 
neutral-current response, because it can be greatly en- 
hanced by coherence, one should extend the study to 
the axial-vector (or spin) response, as it dominates the 
scattering when it is coming from uncorrelated nucleons. 
Further, one should also calculate charged-current inter- 
actions in nonuniform matter. Finally, in the present 
contribution no effort was made to optimize the NN in- 
teraction. While it may be advantageous to do so, any 
accurately calibrated interaction must retain the essen- 
tial features of frustration. Moreover, more sophisticated 
interactions that include momentum and/or density de- 
pendence will significantly increase the computational de- 
mands. At present, we are checking our results against 
more sophisticated simulations using molecular dynam- 
ics, studying finite size effects in larger simulations, ex- 
ploring the temperature and density dependence of our 
results, and calculating the dynamical response. These 
results will be presented in a future contribution |34|. 
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